#!/usr/bin/env python
# -*- coding: UTF-8 -*-
'''
@Project ：V2 
@File    ：dynamic_fjsp_ga.py
@IDE     ：PyCharm 
@Author  ：郭星
@Date    ：2025/9/20 20:30 
'''
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
当然可以——你要的是动态（Dynamic）场景：作业分批到达、机器随机故障/修复，系统滚动重调度。下面给你一份单文件、可直接运行的 Python 代码，包含：

动态 FJSP（柔性作业车间）问题建模

事件：作业释放时间（release）、机器故障（MTBF）/修复（MTTR）（以停机区间的形式）

滚动重调度：每次有事件或工序完成时，用一个小规模 GA对“剩余工序”做全局排程，仅派发最先的一步，然后推进时间（典型 RH/Model-Predictive 控制）

考虑停机区间：解码时遇到停机，自动推迟开始或暂停延长

输出：最终 Cmax、甘特图（带故障窗）、简单日志
dynamic_fjsp_ga.py
---------------------------------
Dynamic Flexible Job Shop Scheduling with Rolling-Horizon GA:
- Job releases (stochastic or preset)
- Machine breakdown/repair windows derived from MTBF/MTTR (downtime intervals)
- Rolling re-optimization: run a small GA on *remaining operations*, dispatch only the first step
- Downtime-aware decoding: task starts avoid downtime; mid-process downtime pauses the task

Outputs:
- Final makespan
- Gantt chart with downtime bands
"""

import math
import random
from dataclasses import dataclass
from typing import Dict, List, Tuple, Optional
import numpy as np
import matplotlib.pyplot as plt

# -------------------------------
# Utilities
# -------------------------------
def seed_all(seed=0):
    random.seed(seed); np.random.seed(seed)

def draw_exponential(mean: float, rng: np.random.RandomState):
    # If mean<=0, return a huge number (no events)
    if mean <= 0: return 1e12
    return float(rng.exponential(mean))

def plot_gantt_with_downtimes(m: int,
                              intervals_by_machine: Dict[int, List[Tuple[float,float,str]]],
                              downtime_by_machine: Dict[int, List[Tuple[float,float]]],
                              title="Dynamic FJSP Gantt"):
    fig, ax = plt.subplots(figsize=(12, 0.8*m + 3))
    yticks, ylabels = [], []
    for mm in range(m):
        # downtime as semi-transparent bars
        for (st, en) in downtime_by_machine.get(mm, []):
            ax.barh(mm, en-st, left=st, color="tab:red", alpha=0.15, edgecolor=None)
        for (st, du, lab) in sorted(intervals_by_machine.get(mm, []), key=lambda x: x[0]):
            ax.barh(mm, du, left=st, color="tab:blue", alpha=0.8)
            ax.text(st + du/2, mm, lab, ha="center", va="center", fontsize=8, color="white")
        yticks.append(mm); ylabels.append(f"M{mm}")
    ax.set_yticks(yticks); ax.set_yticklabels(ylabels)
    ax.set_xlabel("Time"); ax.set_title(title)
    ax.grid(axis="x", alpha=0.3)
    plt.tight_layout(); plt.show()

# -------------------------------
# Problem modeling
# -------------------------------

@dataclass
class Operation:
    machine_time: Dict[int, int]  # feasible machine -> processing time

@dataclass
class Job:
    ops: List[Operation]
    release: float                 # release time of job (time when op0 becomes ready)

@dataclass
class FJSP:
    jobs: List[Job]
    m: int  # number of machines

    @property
    def n_jobs(self): return len(self.jobs)
    @property
    def n_ops_total(self): return sum(len(j.ops) for j in self.jobs)

    def feasible_machines(self, j, k) -> List[int]:
        return list(self.jobs[j].ops[k].machine_time.keys())

# -------------------------------
# Dynamic instance generator
# -------------------------------

def build_dynamic_fjsp(n_jobs=16, m=6, ops_per_job=4, k_machines_per_op=2,
                       pt_low=1, pt_high=9, release_span=40.0,
                       mtbf=120.0, mttr=20.0,
                       horizon=400.0, seed=0) -> Tuple[FJSP, Dict[int, List[Tuple[float,float]]]]:
    """
    - Job releases uniform in [0, release_span]
    - Machine downtimes are generated as alternating up/down periods via Exp(MTBF)/Exp(MTTR)
    """
    rng = np.random.RandomState(seed)

    # Jobs
    jobs=[]
    for _ in range(n_jobs):
        ops=[]
        feas_all = list(range(m))
        for _o in range(ops_per_job):
            feas = sorted(rng.choice(feas_all, size=k_machines_per_op, replace=False))
            mt = {mm: int(rng.randint(pt_low, pt_high+1)) for mm in feas}
            ops.append(Operation(machine_time=mt))
        release = float(rng.uniform(0.0, release_span))
        jobs.append(Job(ops=ops, release=release))
    pb = FJSP(jobs=jobs, m=m)

    # Downtime windows per machine
    downtime = {mm: [] for mm in range(m)}
    for mm in range(m):
        t = 0.0
        while t < horizon:
            up = draw_exponential(mtbf, rng)
            t += up
            if t >= horizon: break
            down = draw_exponential(mttr, rng)
            st = t; en = min(horizon, t + down)
            downtime[mm].append((st, en))
            t = en
    return pb, downtime

# -------------------------------
# Downtime-aware scheduling helpers
# -------------------------------

def adjust_for_downtime(mm:int, start:float, pt:float, downtime_windows: List[Tuple[float,float]]) -> Tuple[float,float]:
    """
    Given a start & processing time on machine mm, return (actual_start, actual_end) accounting downtimes:
    - If start falls in a downtime window => shift to window.end
    - If processing spans into a downtime => pause and resume after it
    """
    t = start
    end_needed = pt
    # if start inside downtime, move to end
    for (st,en) in downtime_windows:
        if t < en and t >= st:
            t = en
            break
    # now accumulate processing while skipping downtimes
    cur = t
    remaining = end_needed
    while remaining > 1e-9:
        # find next downtime that intersects [cur, cur+remaining]
        next_block = None
        for (st,en) in downtime_windows:
            if en <= cur:        # past
                continue
            if st >= cur + remaining:  # future after finish
                break
            if st < cur + remaining and en > cur:  # overlaps
                next_block = (st,en)
                break
        if next_block is None:
            # no overlap -> finish
            return (t, cur + remaining)
        st,en = next_block
        if st <= cur:
            cur = en        # we are in downtime now -> wait
            continue
        # we can process until st
        work = st - cur
        if work >= remaining:
            return (t, cur + remaining)
        else:
            remaining -= work
            cur = en        # resume after downtime
    return (t, cur)

# -------------------------------
# Decoder & GA over remaining ops
# -------------------------------

class RemainingDecoder:
    """
    Decode a plan (op_seq + mac_seq) over remaining operations from the CURRENT state:
    - job_next[j], job_ready[j], mach_ready[mm], downtime
    """
    def __init__(self, pb: FJSP,
                 job_next: np.ndarray, job_ready: np.ndarray,
                 mach_ready: np.ndarray,
                 downtime: Dict[int, List[Tuple[float,float]]]):
        self.pb = pb
        self.job_next = job_next.copy()
        self.job_ready = job_ready.copy()
        self.mach_ready = mach_ready.copy()
        self.downtime = downtime
        self.assignments: Dict[Tuple[int,int], Tuple[int,float,float]] = {}  # for debug/visual

    def decode(self, op_seq: List[int], mac_seq: List[int]) -> Tuple[float, List[Tuple[int,int,int,float,float]]]:
        """
        Returns:
         - makespan
         - a list of scheduled items: (j,k,mm,st,en) in chronological order
        """
        events = []
        job_next = self.job_next.copy()
        job_ready = self.job_ready.copy()
        mach_ready = self.mach_ready.copy()

        mac_ptr = 0
        for j in op_seq:
            k = job_next[j]
            if k >= len(self.pb.jobs[j].ops):
                continue
            feas = self.pb.feasible_machines(j, k)
            if len(feas) == 0:
                continue
            idx = mac_seq[mac_ptr] % len(feas)
            mm = feas[idx]
            pt = float(self.pb.jobs[j].ops[k].machine_time[mm])

            earliest = max(job_ready[j], mach_ready[mm])
            st_eff, en_eff = adjust_for_downtime(mm, earliest, pt, self.downtime.get(mm, []))
            job_ready[j] = en_eff; mach_ready[mm] = en_eff
            job_next[j] += 1
            events.append((j, k, mm, st_eff, en_eff))
            mac_ptr += 1

        mk = float(np.max(job_ready))
        events.sort(key=lambda x: x[3])
        return mk, events

class GA_Remaining:
    def __init__(self, pb:FJSP, decoder_factory, pop=60, gens=120, pc=0.9, pm=0.2, seed=0):
        self.pb=pb; self.decoder_factory=decoder_factory
        self.pop=pop; self.gens=gens; self.pc=pc; self.pm=pm
        seed_all(seed)

    def _base_op_seq(self, job_next: np.ndarray) -> List[int]:
        seq=[]
        for j in range(self.pb.n_jobs):
            rem = len(self.pb.jobs[j].ops) - job_next[j]
            if rem>0: seq += [j]*rem
        return seq

    def rand_ind(self, job_next: np.ndarray) -> Tuple[List[int], List[int]]:
        base = self._base_op_seq(job_next)
        random.shuffle(base)
        mac=[]
        seen=[0]*self.pb.n_jobs
        for j in base:
            k = seen[j]
            feas = self.pb.feasible_machines(j, k)
            mac.append(random.randrange(max(1,len(feas))))
            seen[j]+=1
        return base, mac

    def crossover_ppx(self, a_seq, b_seq):
        n=len(a_seq); need={}
        for x in a_seq: need[x]=need.get(x,0)+1
        left=dict(need); res=[]
        ai, bi = 0, 0
        while len(res)<n:
            pickA = (random.random()<0.5)
            if pickA:
                while ai<n and left.get(a_seq[ai],0)==0: ai+=1
                if ai<n:
                    j=a_seq[ai]; res.append(j); left[j]-=1; ai+=1
                else:
                    while bi<n and left.get(b_seq[bi],0)==0: bi+=1
                    j=b_seq[bi]; res.append(j); left[j]-=1; bi+=1
            else:
                while bi<n and left.get(b_seq[bi],0)==0: bi+=1
                if bi<n:
                    j=b_seq[bi]; res.append(j); left[j]-=1; bi+=1
                else:
                    while ai<n and left.get(a_seq[ai],0)==0: ai+=1
                    j=a_seq[ai]; res.append(j); left[j]-=1; ai+=1
        return res

    def crossover_uniform(self, a, b):
        return [ (aa if random.random()<0.5 else bb) for aa,bb in zip(a,b) ]

    def mutate(self, op_seq, mac_seq, job_next):
        if random.random()<0.5 and len(op_seq)>=2:
            i,j = random.sample(range(len(op_seq)),2)
            op_seq[i], op_seq[j] = op_seq[j], op_seq[i]
        if random.random()<0.7:
            seen=[0]*self.pb.n_jobs
            for i,j in enumerate(op_seq):
                k=seen[j]
                if random.random()<0.2:
                    feas = self.pb.feasible_machines(j, k)
                    mac_seq[i] = random.randrange(max(1,len(feas)))
                seen[j]+=1

    def evaluate(self, op_seq, mac_seq, job_next, job_ready, mach_ready, downtime):
        dec = self.decoder_factory(job_next, job_ready, mach_ready, downtime)
        mk, _ = dec.decode(op_seq, mac_seq)
        return mk

    def run(self, job_next, job_ready, mach_ready, downtime, gens=None):
        gens = gens or self.gens
        pop=[]; fit=[]
        for _ in range(self.pop):
            ind = self.rand_ind(job_next)
            pop.append(ind)
            fit.append(self.evaluate(ind[0], ind[1], job_next, job_ready, mach_ready, downtime))
        for g in range(gens):
            new=[]
            while len(new) < self.pop:
                i,j = random.randrange(self.pop), random.randrange(self.pop)
                p1 = pop[i] if fit[i]<fit[j] else pop[j]
                i,j = random.randrange(self.pop), random.randrange(self.pop)
                p2 = pop[i] if fit[i]<fit[j] else pop[j]
                c1=(p1[0][:], p1[1][:]); c2=(p2[0][:], p2[1][:])
                if random.random()<self.pc:
                    c1=(self.crossover_ppx(p1[0],p2[0]),
                        self.crossover_uniform(p1[1],p2[1]))
                    c2=(self.crossover_ppx(p2[0],p1[0]),
                        self.crossover_uniform(p2[1],p1[1]))
                if random.random()<self.pm: self.mutate(*c1, job_next)
                if random.random()<self.pm: self.mutate(*c2, job_next)
                new += [c1,c2]
            pop = new[:self.pop]
            fit = [self.evaluate(ind[0], ind[1], job_next, job_ready, mach_ready, downtime) for ind in pop]
        idx = int(np.argmin(fit))
        return pop[idx], fit[idx]

# -------------------------------
# Rolling-horizon simulator
# -------------------------------

class DynamicSimulator:
    def __init__(self, pb:FJSP, downtime:Dict[int,List[Tuple[float,float]]], seed=0):
        self.pb = pb
        self.downtime = downtime
        self.m = pb.m
        seed_all(seed)

        self.time = 0.0
        self.job_next = np.zeros(pb.n_jobs, dtype=int)
        self.job_ready = np.array([job.release if len(job.ops)>0 else 1e18 for job in pb.jobs], dtype=float)
        self.mach_ready = np.zeros(self.m, dtype=float)

        self.done_ops = 0
        self.total_ops = pb.n_ops_total

        self.timeline: Dict[int, List[Tuple[float,float,str]]] = {mm:[] for mm in range(self.m)}  # per machine: (start, dur, label)

    def remaining_exists(self):
        return self.done_ops < self.total_ops

    def next_external_event_time(self):
        # The next job release that is > current time and whose job still has ops remaining
        tmin = 1e18
        for j, job in enumerate(self.pb.jobs):
            if self.job_next[j] < len(job.ops):
                r = job.release
                if r > self.time:
                    tmin = min(tmin, r)
        return tmin

    def all_ready_jobs(self):
        ready = []
        for j in range(self.pb.n_jobs):
            if self.job_next[j] < len(self.pb.jobs[j].ops):
                if self.job_ready[j] <= self.time + 1e-9:
                    ready.append(j)
        return ready

    def step(self, ga_pop=50, ga_gens=80):
        """
        One RH step:
         - If nothing is ready, jump to next release time (or earliest machine ready if both future).
         - Else run GA on remaining ops, decode to get a full plan; dispatch only the FIRST finishing operation from that plan.
        """
        # If no job is ready now, fast-forward
        ready_jobs = self.all_ready_jobs()
        if len(ready_jobs) == 0:
            t_rel = self.next_external_event_time()
            t_mr  = float(np.min(self.mach_ready))
            t_next = min(t_rel, t_mr) if t_rel < 1e17 else t_mr
            self.time = max(self.time, t_next)
            return

        # Build GA on remaining ops
        decoder_factory = lambda jn, jr, mr, dt: RemainingDecoder(self.pb, jn, jr, mr, dt)
        ga = GA_Remaining(self.pb, decoder_factory, pop=ga_pop, gens=ga_gens, seed=random.randrange(1<<30))
        (op_seq, mac_seq), best_val = ga.run(self.job_next, self.job_ready, self.mach_ready, self.downtime, gens=ga_gens)

        # Decode once to get the plan & events
        dec = decoder_factory(self.job_next, self.job_ready, self.mach_ready, self.downtime)
        _, plan = dec.decode(op_seq, mac_seq)
        if len(plan) == 0:
            # Nothing can be scheduled right now, advance time a little (to next release/machine)
            self.time = max(self.time, min(self.next_external_event_time(), float(np.min(self.mach_ready))))
            return

        # Pick the *soonest finishing* operation from plan that *starts not later than now+eps* if possible,
        # else take the earliest-start operation.
        eps = 1e-6
        feasible_now = [(j,k,mm,st,en) for (j,k,mm,st,en) in plan if st <= self.time + eps]
        if len(feasible_now) == 0:
            # Take the earliest start item; fast-forward time to its start
            (j,k,mm,st,en) = sorted(plan, key=lambda x: x[3])[0]
            self.time = max(self.time, st)

        # Recompute feasible-now with updated time
        feasible_now = [(j,k,mm,st,en) for (j,k,mm,st,en) in plan if st <= self.time + eps]
        if len(feasible_now) == 0:
            # Still none; small nudge
            self.time += 1e-3
            return

        # Execute the one with earliest *end*
        (j,k,mm,st,en) = sorted(feasible_now, key=lambda x: x[4])[0]
        duration = en - st
        # Record on timeline
        self.timeline[mm].append((st, duration, f"J{j}-O{k}"))
        # Advance state
        self.time = en
        self.job_next[j] += 1
        self.job_ready[j] = en
        self.mach_ready[mm] = en
        self.done_ops += 1

    def run(self, verbose=True, max_steps=100000):
        steps=0
        while self.remaining_exists() and steps < max_steps:
            self.step()
            steps += 1
            if verbose and steps % 100 == 0:
                done_pct = 100.0 * self.done_ops / self.total_ops
                print(f"[Sim] time={self.time:.2f}  done={self.done_ops}/{self.total_ops} ({done_pct:.1f}%)")
        makespan = self.time
        return makespan, self.timeline

# -------------------------------
# Demo
# -------------------------------

def main():
    seed_all(0)
    # Build a dynamic FJSP instance
    pb, downtime = build_dynamic_fjsp(
        n_jobs=18, m=6, ops_per_job=4, k_machines_per_op=2,
        pt_low=1, pt_high=9,
        release_span=40.0,
        mtbf=120.0, mttr=20.0,
        horizon=400.0,
        seed=7
    )

    sim = DynamicSimulator(pb, downtime, seed=1)
    mk, timeline = sim.run(verbose=True)
    print(f"\nFinal makespan (dynamic) = {mk:.2f}")

    # Plot Gantt with downtime bands
    plot_gantt_with_downtimes(pb.m, timeline, downtime, title=f"Dynamic FJSP (Cmax={mk:.1f})")

if __name__ == "__main__":
    main()
